% Computes bi-variate transfer functions between two channels, given
% two input & one output channel
% Latest date 21.11.2002
% bug fixed - correct signs for the cross spectra
% compared with readedi1 erros, z coh found correct
% introduced flexible array (for input) 22.3.2003

%zij = za zb - zc zd;
%tf = zxy zyx zxx zyy

function[t] = tf1(SPMdata,ProcDef),

A = size(SPMdata);
len = length(A);
p = sqrt(-1);
nfreq=ProcDef.nfrq;
EQN=ProcDef.equations;
neq=ProcDef.neq;


if len == 2,
   SPM(1,:,:) = SPMdata(:,:);
   t = zeros([1 neq]);
elseif len == 3,
   SPM(:,:,:) = SPMdata(:,:,:);
   t = zeros([A(1) neq]);
elseif len == 4 & A(1) == 1 & A(2) == 1,
   SPM(1,:,:) = SPMdata(1,1,:,:);
   t = zeros([1 neq]);
elseif len == 4 & A(1) == 1,
   SPM(:,:,:) = SPMdata(1,:,:,:);
   t = zeros([A(2) neq]);
elseif len == 4 & A(2) == 1,
	SPM(:,:,:) = SPMdata(:,1,:,:);   
   t = zeros([A(1) neq]);
else,
   t = 0;
   fprintf('tf.m -> Cannot handle this arrays\n');
   return;
end;


for l = 1:neq,
i = EQN(l,1);
j = EQN(l,2);
k = EQN(l,3);
za = SPM(:,j,j);
if i>k,
   zb = SPM(:,i,k)+p*SPM(:,k,i);
elseif k>i,
   zb = SPM(:,k,i)-p*SPM(:,i,k);
end;

if k>j,
   zc = SPM(:,k,j)-p*SPM(:,j,k);
elseif k<j,
   zc = SPM(:,j,k)+p*SPM(:,k,j);
end;

if i>j,
   zd = SPM(:,i,j)+p*SPM(:,j,i);
elseif i<j,
   zd = SPM(:,j,i)-p*SPM(:,i,j);
end;

den = SPM(:,j,j).*SPM(:,k,k)- (SPM(:,k,j) - p*SPM(:,j,k)).*(SPM(:,k,j) + p*SPM(:,j,k));
top = [za.*zb-zc.*zd SPM(:,k,k).*zd-conj(zc).*zb];

t(:,l)   = (top(:,1)./den);
t(:,l+2) = (top(:,2)./den);
end;
